Bose-glass to Superfluid transition in the three-dimensional Bose-Hubbard Model 
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We present a Monte Carlo study of the Bose-glass to superfluid transition in the three-dimensional 
Bose-Hubbard model. Simulations are performed on the classical (3 -f- 1) dimensional link-current 
representation using the geometrical worm algorithm. Finite-size scaling analysis (on lattices as large 
as 16x16x16x512 sites) of the superfluid stiffness and the compressibility is consistent with a value 
of the dynamical critical exponent 2; = 3, in agreement with existing scaling and renormalization 
group arguments that z = d. We find also a value of 1/ = 0.70(12) for the correlation length 
exponent, satisfying the relation v > 2/d. However, a detailed study of the correlation functions, 
C(r, r), at the quantum critical point are not consistent with this value of z. We speculate that 
this discrepancy could be due to the fact that the correlation functions have not reached their true 
asymptotic behavior because of the relatively small spatial extent of the lattices used in the present 
study. 

PACS numbers: 67.40.Db, 67.90.+Z 



I. INTRODUCTION 

Quantum phase transitions as they occur in models 
comprised of bosons have been the focus of considerable 
interest lately. Most notably, experiments by Greiner et 
al.^ have observed the transition between a Mott insulat- 
ing (MI) phase and a superfluid phase in an optical lat- 
tice loaded with ^^Rb atoms. The observed phase tran- 
sition between the superfluid and the insulating phase 
is thought to share the universal properties of a variety 
of physical systems, including He'' in porous substrates^, 
Josephson-junction arrays^ as well as thin superconduct- 
ing films^i^. While disorder can be neglected in the 
experiments by Greiner et alX, it clearly plays a central 
role in several of the experiments just mentioned. As out- 
lined in the seminal paper by Fisher et alJ the statics as 
well as the dynamics of the quantum phase transitions 
occurring in the Bose-Hubbard model are strongly influ- 
enced by disorder and a new insulating glassy phase, the 
Bose-glass phase, should occur. The transition of inter- 
est in this case is directly between the superfluid (SF) 
and the Bose-glass (BG). Recent experiments showing 
that this disordered transition can also be studied using 
optical latticed has therefore created considerable ex- 
citement. Previous theoretical studieaSiiSiiiiiSiiii^ have 
mainly focused on the transition as it occurs in one or 
two dimensions and very recent theoretical worksi^ti^, 
addressing directly the experimental situation relevant 
for optical lattices, have also mainly analyzed this case. 
In light of the experiments by Lye et al&i^, we consider 
in the present study the transition as it occurs in three 
dimensions in the presence of short-range interactions. In 
particular we determine the dynamical critical exponent, 
z, characterizing this transition in d = 3. 

One hall-mark of the Mott insulating phase is that it is 
incompressible. In contrast to this, both the superfluid 
and the Bose-glass phase are compressible and have a 
flnite non-zero compressibility, k ^ 0. Close to the quan- 
tum critical point between the Bose-glass and the super- 



fluid scaling argumentsi, based on generalized Josephson 
relations, then show that: 



(1) 



where S describes the distance to the critical point in 
terms of the control parameter. Since both the BG and 
SF phases are compressible it then seems natural to ex- 
pect the compressibility to be flnite also at the critical 
point and one must then conclude that 



= d. 



(2) 



It can in fact be shown that assuming the compressibil- 
ity cither diverges or tends to zero at the BG-SF critical 
point leads to implausible behavior. In Ref.0the relation 
z — d was therefore argued to hold not only below the 
upper critical dimension but for any dimension d > 1. 
This result implies that there is not an upper critical di- 
mension for this transition in any conventional sense. In 
d = 1 analytical workS strongly supports the conclusion 
that z = d = 1 and numerical workiiii^iiS in d = 2 also 
flnd z — d = 2. Long-range interactions are likely to 
yield a different ai, but we shall not be concerned with 
that case here. Based on Dorogovtsev'siS double-e ex- 
pansion, an alternative scenario has been proposed^SiSi 
in which the critical exponents jump discontinuously to 
their mean field values at the critical dimension dc — 4. 
In this scenario, for d > dc two stable fixed points are 
present: the gaussian fixed point is stable at weak dis- 
order with the random fixed point remaining stable at 
strong disorder. Below the critical dimension, d < dc, 
only the random flxed point is stable. It should be noted 
that in order to obtain sensible answers from the double-e 
expansion an ultraviolet frequency cutoff uj\ has to be in- 
troducedSS — a procedure which seems difficult to justify. 
This was later remedied uponSi still within the frame- 
work of the double-e expansion. It is also possible to 
question the validity of any e-expansion approach on the 
grounds that the existence of an upper critical dimension 



2 



for the disordered transition is implicitly assumed. One 
might then ask if bounds on the upper critical dimensions 
exist. Indeed, using the exact inequalitjiS^ 



>2M 



(3) 



valid for stable fixed points in the presence of correlated 
disorder, one immediately sees that the requirement that 
u = 1/2 at the upper critical dimension, dc, yields dc > 4. 
Since the existence of an upper critical dimension is de- 
batable it is then natural to instead focus on the lower- 
critical dimension, which is well established. Indeed, it 
is also possible to develop an expansion away from the 
lower-critical dimension {di = Using this ap- 

proach the dynamical critical exponent is exactly z — d 
and the correlation length exponent v can be calculated 
to second order in \/c? — 1 yielding good agreement with 
experimental and numerical results in d = 2. The onset 
of mean field behavior, if any, is therefore at best uncon- 
ventional in this model and still a matter of debate. In 
particular, it would be valuable to know if the relation 
z — d continues to hold in dimensions higher than d = 2. 

In the present work, we present a Monte Carlo study 
of the three-dimensional Bose-glass to superfluid transi- 
tion at strong disorder. Our focus is reliable estimates 
of the dynamical critical exponent z in order to test 
the relation z = d = 3. As outlined above, much of 
the numerical work to date has focused on the one- and 
two-dimensional cases, which are less demanding compu- 
tationally. However, recently a very efficient geometric 
worm algorithmi2i2^ has been developed for the study of 
bosonic phase transitions making significantly larger sys- 
tem sizes available. Even though the geometrical worm 
algorithm we use has proven to be very efficient for deal- 
ing with large lattices in two dimensions (also in the pres- 
ence of disorder), we were only able to study a relatively 
limited number of system sizes in d = 3 that were large 
enough to avoid finite-size effects but small enough to 
properly equilibrate with a feasible amount of computa- 
tional effort. In the present study therefore, we focus 
exclusively on the three dimensional case. We leave for 
future study the case of d = 4 that would be of consider- 
able interest for the approach to mean field suggested by 
Weichman and collaborators22i2i. 

It has been suggested that in the vicinity of commen- 
surate values of the chemical potential disorder is weakly 
irrelevant2Si2L2Si2S. In this case a direct MI-SF transition 
could occur even in the presence of weak disorder. How- 
ever, recent high precision numerical work^ii reached the 
opposite conclusion for d = 2, although very large sys- 
tem sizes were needed in order to show this. Since we are 
mainly interested in the BG-SF transition, we minimize 
any cross-over effects by performing all of our calcula- 
tions at a value of the chemical potential where in the 
absence of any disorder the system is in the superfluid 
phase for any non-zero hopping and there is no MI-SF 
transition^. The BG-SF transition we observe is there- 
fore induced by the disorder and corresponds in the RG 
sense to a random fixed point. 



The paper is outlined as follows: the Bose-Hubbard 
model is introduced below in further detail, including the 
mapping to a (d-|-l) dimensional classical model on which 
we perform our study. Section IlII Al outlines the scaling 
theories necessary to extract the critical exponents. Sec- 
tion ^^B] details the numerical techniques and discusses 
the difficulties we encountered in properly equilibrating 
our simulations. Finally, we present our results in section 

El 



II. MODEL 

We begin with the Bose-Hubbard Hamiltonian, includ- 
ing on-site disorder in the chemical potential»i 



E 



{r,r') 



(*J;*r'+H.C. 



(4) 

Here and $r are boson creation and annihilation op- 
erators at site r, and = ^J^r is the number operator. 
The on-site repulsive interaction U localizes the bosons 
and competes with the delocalizing effects of the tunnel- 
ing coefficient t. The random chemical potential fir is 
distributed uniformly on [/x — A, /u -I- A]. As noted by 
Damski et alm^, the introduction of a random potential 
in an optical lattice will also generate randomness in the 
hopping (tunneling) term, t. However, this type of disor- 
der can be ignored^®. The phase diagrami for the pure 
model consists of a superfluid phase at high t/U which is 
unstable to a series of Mott-insulating regimes centered 
at commensurate densities at low t/U. In the presence 
of disorder, the Bose-glass phase stabilizes between the 
insulating and superfluid phases. 

Following standard methodsii*^, by integrating out 
amplitude fluctuations of the Bose field to second order, 
we transform the Bose-Hubbard Hamiltonian into an ef- 
fective classical Hamiltonian in {d + 1) dimensions well 
suited for Monte Carlo study: 



H 



K ^ 12 



(r,r) 



(5) 



(r,r) 



The integer currents J(r,T) ^re defined on the bonds of the 
lattice and obey a divergenceless constraint VJ(r_T-) — 0. 
The resulting current loops are interpreted in this context 
as world lines of boson^ii*^ and represent fluctuations 
from an average, non-zero density. The coupling K acts 
as an effective classical temperature and drives the phase 
transition: at low K the repulsive short range interaction 
U dominates and the system is insulating, while at high 
K the hopping t dominates and the system is superfluid. 
In this transformation amplitude fluctuations of the bo- 
son fields are integrated out to second order. However, 
these fluctuations are not expected to affect the universal 
details of the transition. 

The advantage of the formulation of the Bose-Hubbard 
model in terms of the link-currents, Eq. |SJ|, is that an 
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extremely ef&cient worm algorithmiSi is available for this 
model. This worm algorithm can also be directed^ even 
in the presence of disorder. However, the memory over- 
head associated with using a directed algorithm is pro- 
hibitive for the present study due to the presence of dis- 
order and the high spatial dimension. We have therefore 
exclusively used the more standard formulation of the 
algorithmic. 

III. METHOD 
A. Scaling 

The BG-SF transition is a continuous quantum phase 
transition, and is characterized by two diverging length 
scales: a spatial (^) and a temporal (fr) correlation 
length, related by the dynamical critical exponent z: 

^ r ~ (5--)-, (6) 

with (5 = [K — Kc)/Kc. These form the basis of the 
scaling theory used to analyze our data. 

The two observables of primary interest are the super- 
fluid stiffness p and the compressibility k. The super- 
fluid stiffness p (proportional to the superfluid density) 
is defined in terms of the change in free energy associ- 
ated with a twist in the spatial boundary conditions. As 
for the compressibility, the critical behaviour of p can be 
derived as a generalization of the Josephson scaling re- 
lations for the classical transitionLiiiSi, and scales with 
the correlation length as 

^_^-(d+.-2)^ (7) 

The compressibility can similarly be found as a re- 
sponse to a twist applied to the temporal boundary con- 
ditions, and is found to scale as^iLSi 

«^r''-^\ (8) 

leading to the relation z = d previously mentioned. 

The first step in the study of the critical properties of 
the Bose- Hubbard model is a precise determination of the 
location of the critical point through finite-size scaling 
analysis. The presence of two correlation lengths implies 
that finite-size scaling functions will have two arguments: 
f{L/^,Lr/^r)- Hence 

p = ^-(d+^-^)f{L/tLr/^r), (9) 

or, by appropriately scaling the arguments and requiring 
that the stiffness remain finite at the critical point for a 
system of finite size, 

P^J^I^^PiL'/'S^Lr/Ln. (10) 

This complicates the scaling analysis as we must work in 
a two-dimensional space. The first approach is to work at 



lattice sizes whose temporal sizes scale with the exponent 
z; that is to work at a fixed aspect ratio: 

a^Lr/L\ (11) 

The quantity: 

L^+^-^piL'/'^S,a^^) (12) 

should then be a universal function of a at Kc- We can 
hold the aspect ratio constant by working with systems 
of dimension L'^ x aL^ . If an initial estimate for the dy- 
namical critical exponent z is available, the critical point 
can be located by plotting L'^^^~'^ p{L^/^ 5, a) versus K 
for several different linear system sizes L. If the correct 
value of z is used, these curves will all intersect at the 
critical point. This unfortunately requires that the initial 
estimate for the value of z be made before simulations are 
run. 

Alternatively, with an estimate of Kc^ the first argu- 
ment of Eq. E|can be held constant. Curves of pL'^^''-~^ 
plotted for different L against the ratio Lr/L^ should 
then collapse for the correct value of We use 

both approaches, with the hope that a consistent picture 
emerges. 

Once Kc has been located by the above method we can 
proceed to study the behavior of the compressibility at 
the critical point. If indeed z is equal to d, as described in 
the introduction, the compressibility should be roughly 
constant as K varies through K^^ and should not show 
any dependence on L. In particular k should neither 
diverge nor go to zero at Kc- 

We also consider the particle-particle correlation func- 
tion C(r,T), which is expectedi to decay asymptotically 
as 

C(r, T = 0) - r-y^ , C(r = 0, r) - r"'^- (13) 
with exponents given by 

2/r = d + Z - 2 + 7], 

Vr = (d + z-2 + r;)/z. (14) 

In addition to defining the exponent 77, with reliable esti- 
mates of the correlation functions this would in principle 
provide for an independent calculation of the dynamical 
critical exponent z. Fisher et al7 also derive bounds for 
77: 

2- {d + z) <r] <2~ d. (15) 

The lower bound arises from the requirement that the 
correlation functions decay, and the upper bound is ar- 
gued from the scaling of the density of states in the Bose- 
glass and superfluid phases. These can be simply stated 
as bounds on the exponents yr and y-r'. 

< yr < 1, Q<yr<z. (16) 
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FIG. 1: Main panel: histograms of {pL"^) for a set of 1000 dis- 
order realizations on a system of size 8x8x8x64 at our estimate 
of Kc = 0.190(1). Curves show the evolution of the distribu- 
tion P{{pL'^)) as a function of the number of Monte Carlo 
sweeps ts performed after equilibration on each disorder real- 
ization. The peak in the distributions at {pL"^) = persists 
for many times the relaxation time, as measured by the con- 
vergence of the Hamming distances, while the broader peak 
near (pL*) = 0.025 grows. Inset: Hamming distances calcu- 
lated on the same set of disorder realizations. Ha.ji is plotted 
against the total number of sweeps performed, t. Ha.ag is 
plotted against the number of sweeps performed after equili- 
bration, ta- Open symbols are calculated in terms of spatial 
currents. Solid symbols are calculated in terms of temporal 
currents. 

The final exponent we calculate is the correlation 
length exponent ly. Taking the derivative of (|I0|) with 
respect to K we see that 

^d+z-2dp_ ^ l^l/^plfj^l/.^^ Lr/L'). (17) 

dK 

Plotting this derivative against L at Kc should yield 
a power law with exponent l/v. The crossing data 
pj^d+z-2 calculated for different system sizes with a fixed 
aspect ratio, a, should also collapse to a single curve 
when plotted against L^^^S. 

B. Numerical Method 

Both the superfluid stiffness and the compressibility 
can be calculated in terms of the link-current wind- 
ing numbersiii^ = L^^J^tt^Zt ea.ch direction 
^ — x,y, z,T (here L^^y^z — L). The superfluid stiffness 
is related to a twist in the spatial boundary conditions 
and so can be calculated in terms of the spatial winding 
numbers: 

P= j^rf\j ("7=x,y,.)],.,- (18) 

We denote thermal averages by (O) and disorder aver- 
ages by [C]^^. Similarly, the compressibility is associ- 
ated with a twist in the temporal boundary conditions 



FIG. 2: Distributions of (pL*) for L = 8 and 12 at A'c, plotted 
on linear (main panel) and logarithmic (inset) axes. The av- 
erage of each distributions is close to [{pL'^)\av = 0.08, which 
is significantly higher than the typical values of the distribu- 
tions, indicating that the broad tail of the distribution must 
be well sampled to obtain an accurate estimate of the true 
average. 

and is defined in terms of the temporal winding num- 
bers: 

^=^[{nl)~{nr)o.{nr)p\^,. (19) 

Since this expression contains the disorder average of the 
square of a thermal average, the systematic error is re- 
duced by calculating this average on two independent 
lattices a and (3 with the same disorder realization^iS^. 

We are also interested in the derivative of p with re- 
spect to the coupling. This can be calculated thermody- 
namically from the stiffness and total energy E: 

^-^.[^PE)-{P)^{E),]^^. (20) 

This was found to produce better estimates than numer- 
ically differentiating p in the two-dimensional case^*i. 

Finally, since the construction of each 'worm' is es- 
sentially equivalent to propagating a boson through the 
lattice, it is possible to calculate the particle-particle cor- 
relation function directly from the behaviour of the algo- 
rithm. Details can be found in Ref. [2^ 

As always, the system must be be run for Monte 
Carlo sweeps at each disorder realization to ensure that 
equilibrium has been reached before beginning to sample 
the generated configurations. To confirm this, two sim- 
ulations are carried out simultaneously on lattices with 
the same disorder realization but different initial config- 
urations (we set all of the currents in each direction to a 
different integer constant for a and (3) . These initial con- 
figurations are far from equilibrium. It is useful to define 
'Hamming distances' between different current configu- 
rations in order to measure the relaxation time of the 
algorithm (this is done in the spirit of Ref. |^ though 
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FIG. 3: The main panel shows pL**^^^^ = pL* plotted versus 
K for lattices of size Lr = assuming d = z = 3. The 
lines are guides to the eye. The crossing gives an estimate of 
the critical point Kc = 0.190(1) and is consistent with z = 
3. The two insets show equivalent results with two different 
aspect ratios a — 1/32 and a = 1/64. 

the definitions used here are shghtly different). We de- 
fine the Hamming distance between the two lattices a 
and f3 after performing a total of t Monte Carlo sweeps 
on their initial configurations: 

^-/"W = ^ [^«,(r,.)W - Jhr.m"- (21) 

Similarly, we define the Hamming distance between the 
configuration of lattice a at the sweep where we begin 
to sample the generated configurations (denoted ao) and 
the configuration of a after a further tg — t ~ to sweeps: 

Ha'ao'^iis) = Jjr[~ X! [•^a,(r,T){io+ts) - -'"^.(r.r) (^o)] • 

In the present study, Iq is chosen prior to running the 
simulations and is held constant. 

For the sake of simplicity, we define a Monte Carlo 
sweep to be the construction of a single worm; while 
not ideal for comparing its characteristics to other al- 
gorithms, this is sufficient for our discussion here. Since 
the initial configuration of a and (3 are different, iJ^ pi^) 
will be large at the beginning of the simulation [t small) . 
Conversely, since shortly after the configuration of a 
will not have changed substantially, H^^^^{ts) will be 
small. If to has been chosen larger than the relaxation 
time of the algorithm, then for sufficiently large values 
of t and ta the configurations of ao, a, and [3 will be 
independent, equilibrated states, and thus the two dis- 
tances should converge. The inset of figure shows 
this approach to equilibrium on an 8x8x8x64 lattice with 
to = 3 X 10^ at our estimate of Kc- The distances con- 
verge after t « 10^, indicating that the relaxation time 
for this system is approximately thirty thousand sweeps. 

A further complication in the calculation of the dis- 
order averages was encountered due to the discrete na- 
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FIG. 4: Search for a crossing in pL'^'^'~^ for z — 2 using 
lattices of size x ai^^'^, (see also Fig. EJ. Results are 
shown for two different aspect ratios a = 1/8, 1/16. The 
crossings show significant drift between different lattice sizes 
and between the two aspect ratios shown. Consequently, z ^ 
2. 

turc of the winding numbers which are used to calculate 
p and K. Since most disorder realizations have a small 
but finite superfluid stiffness at Kc, many independent 
configurations of must be generated for each disorder 
realization to achieve a reliable estimate. Figure^shows 
the distribution of {pL^) as a function of the number of 
sweeps ts performed after equilibration on each realiza- 
tion of the disorder. For t^ « t^, there are many realiza- 
tions for which no configuration is generated with a finite 
winding number. Only after running for much longer at 
each realization are the true features of the distribution 
resolved. Under-sampling these realizations runs the risk 
of underestimating the disorder average. This issue is 
discussed at further length in Ref. 

The computational demands of equilibration grow 
quickly with system size. Typical systems we studied, 
of linear size L = 10 to L = 14, required up to 3 x 10^ 
sweeps to relax and up to a further 1 x 10^ sweeps for the 
distributions of {pL'^) to equilibrate. For L — 16, nearly 
1 X 10^ sweeps were required to relax and we were unable 
to run simulations for sufficiently long to see the distri- 
butions equilibrate. Some results for L ~ IQ are shown 
below as they demonstrate consistent scaling. Figure [21 
shows the equilibrated distributions of {pL^) for two lat- 
tice sizes, L = 8 and 12. They show the broad tail which 
necessitates running on at least 10'^ disorder realizations. 
Moreover, at least for the system sizes we were able to 
study, the distributions do not narrow as L increases. 
Hence, the system is likely not self- averaging^ 



IV. RESULTS 

We begin by a discussion of our results for the scal- 
ing of the stiffness p. All our results are for the three- 
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FIG. 5: pL* versus a = Lr/L'' for system sizes L = 8, 10, 12 
assuming z = 3 at if = 0.19. Lattices of size x Lt 
were used. The inflection points and curvature show a very 
good coUapse indicating that the dynamical critical exponent 
is likely 2 = 3. The slight vertical drift with increasing L 
likely indicates that we are slightly off the true Kc = 0.190(1), 
within the bounds of our errorbars. 



dimensional case, hence, the scaling relation Ea.ll2lstates 
that L^^^p{L^/''5, a) should display a crossing at the crit- 
ical point. As already mentioned, we focus on the value of 
the chemical potential /ir distributed randomly between 

— A, /i + A] with and A = i. 

We begin with the ansatz z = d = 3, using lattices of 
size L X L X L X aL^^^. In Fig. |31we show plots of pL"^ 
versus K for linear system sizes ranging from L = 8 to 
16. A clean crossing is observed at Kc = 0.190(1) for 
three values of the aspect ratio, a = 1/8 (main panel), 
a = 1/32, and a = 1/64 (insets). As can be seen from 
Fig. 13 our results for systems of size L = 8 do not scale 
as well as larger system sizes. This effect is more pro- 
nounced for L < 8 (not shown). This is most likely due 
to corrections to scaling which for these small system 
sizes can not be neglected. We also note that an im- 
provement in the scaling behaviour of L = 8 with the 
aspect ratio is visible in Fig. |3| This is presumably due 
to the fact that the optimal value for a is given by the 
relation (^/L)^ ~ ^r/L^. 

In order to test how sensitive we are to the value of z 
we also tried z = 2, the value for the dynamical critical 
exponent quite well established in two dimensions. Our 
results for this value of z are shown in Fig. ^ for lattices 
of dimension x aL^^^. For this value of z our results 
show significant drift in the crossing for different system 
sizes and different aspect ratios as can be clearly seen in 
Fig. 2] The apparent crossing between two system sizes 
are only slightly higher than the clear crossing seen using 
z = 3 in Fig. 13 but only two system sizes can be made to 
cross at a given K. From the results in Fig.^we conclude 
that z^2. 

As explained in section FlII Al we can now hold constant 
the first argument of the scaling function (|12|l by working 
at our estimate of Kc- Assuming that our estimate of 



FIG. 6: The compressibility, k shown versus K near Kc = 
0.190(1). Resuhs are shown for L = 8, 10, 12, and 14 with 
a — 1/8. Note the absence of features at Kc and the lack of 
dependence on system size, indicating again that z = d — 3. 



Kc — 0.190(1) obtained with z = 3 is correct we can 
now check these estimates self-consistently by plotting 
pj^d+z-2=A ygj-g^s lj^jlj'-='i for a range of L at Kc — 
0.190(1). For a detailed explanation of the procedure we 
refer to Ref [S^. Our results are shown in Fig. The 
curves show /5L'^+^~^='* plotted for systems of linear size 
L = 8, 10, and 12 as a function of the aspect ratio a — 
Lr/L^^^- All the curves for different L and Lr collapse 
onto a single curve for 2 = 3. This result is a rather 
strong confirmation that the assumption z = 3 indeed is 
correct. If one studies the curves in detail a very slight 
vertical drift with increasing L is noticeable. This drift is 
likely the result of a small deviation, within our errorbars, 
of the actual critical temperature from our estimate of 
Kc = 0.190(1). It would have been quite interesting to 
study systems of linear size L > 12 or systems with an 
aspect ratio largely exceeding Lr/L^ ~ 1, however, given 
our value of the dynamical critical exponent of z = 3 this 
is computationally too demanding. 

We now discuss our results for the compressibility, k. 
In Fig. the compressibility is shown for a range of K 
around the estimated critical point Kc = 0.190(1) for 
several different system sizes, using the same lattice sizes 
as for in Fig. |31 As expected if d = z, the compressibility 
shows no dependence on the linear size of the system and 
is approximately constant through the transition. Fur- 
thermore, at Kc the compressibility neither diverges nor 
does it tend to zero. This is again a strong confirmation 
that z = d = 3. One should note that the results in 
Fig. are only really useful for determining z if the criti- 
cal point, Kc is known (since one generally would expect 
K to be independent of L far away from the critical point 
where £, ^ L). In the absence of any knowledge of Kc, k 
would have to be calculated for the entire range of phys- 
ically relevant K. It is therefore interesting to compare 
the results shown in Fig. El with the attempt of locating 
the critical point assuming z = 2 shown in Fig. ^ The 
latter results show curves intersecting pairwise around 
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FIG. 7: The derivative of the stiffness, L*dp/dK, as calcu- 
lated directly during the simulation, plotted against the lin- 
ear system size L. Data is shown for a range of couplings 
near Kc- The solid lines indicate power-law fits to the data 
a,t K = 0.190, with exponent 1.11 and K = 0.191 with an 
exponent of 1.69. Data is also shown for simulations run in 
the micro-canonical ensemble of disorder at K — 0.190 (see 
text). The dashed line indicates a power-law fit to the micro- 
canonical data with an exponent indistinguishable from the fit 
to the canonical data. The inset shows the best collapse of the 
scaling curves pL* plotted against SL^^", with Kc = 0.1902 
and u = 0.70. Combining these results we find v = 0.70(12). 
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FIG. 8: Spatial and temporal dependence of the particle- 
particle correlation function, C(r,r — 0) and C(r — 0,r). 
Calculated at K ^ 0.19 and a = 1/8 for L = 8 to 16. Solid 
lines are a fit to C(r,T = 0) = A{r''"'' + {L - r)"^--) and 
C(r = 0,t) = A{t~^^ -I- {L^ - ry^^). Each correlation func- 
tion was fitted with the same power, but a unique coefficient. 
Fitted values are yr — 4.4(4) and j/^ = 1.06(7). Also shown 
are curves calculated for a system of size 20^ x40(a = 1/200) 
at K = 0.1905. This correlation function shows a clear expo- 
nential dependence in the spatial direction, though the tem- 
poral dependence is consistent with the same power law decay 
as for systems of larger aspect ratio. 



K ^ 0.2 with the intersections moving downwards with 
L. If indeed the dynamical critical exponent was z = 2 
and not 2: = 3 as we show here, then it would seem very 
unlikely that the compressibility shown in Fig. (SJ just 
below this range of K ^ 0.2, could be so featureless. 
One would in that case have expected it to behave as 
K ^ S^id.''^)='^ ^ resulting in the finite-size scaling form 
K = (1/L''-^=1)k(L/^, Lr/L''). Since the results in Fig.El 
are independent of L they therefore clearly exclude z = 2. 

We now turn to our results for the correlation length 
exponent v. Figure [7| shows in the main panel a plot of 
the derivative of the stiffness (times L^) for the system 
sizes we studied. This quantity is calculated during the 
simulation without the use of any numerical derivatives. 
As explained in section UlI Al we expect this derivative to 
behave as - L^/" at K^. The results for K = 0.188 
clearly deviate from a power-law consistent with this 
value of K being below Kc- From the results shown in 
Fig. Owe know that K = 0.191 is likely above Kc, while 
K = 0.190 is likely very close to Kc- Hence we fit the 
results for K — 0.190 and K — 0.191 to a power-law, 
L^dp/dK = cL^^'^ , finding an exponent of 1.11 and 1.69, 
respectively (shown as the solid lines in Fig. [TJ. This 
allows us to bracket the estimate of ly to the interval 
0.59 — 0.91. The rather broad range of this interval is 
due to the fact that the relatively large value of z makes 
this way of determining extremely sensitive to a precise 
determination of Kc- We therefore combine this estimate 
with a standard scaling plot of pL versus SL^/" shown 
in the inset of Fig. [7| The best scaling plot is obtained 
with 1^ = 0.70 and Kc = 0.1902 (well within the er- 



rorbars of our estimate for Kc). Combining these two 
estimates of 1/ we conclude v = 0.70(12). This value of 1/ 
satisfies the "quantum" Harris criterion, v > 2/d, which 
should hold for all disorder-driven transitions^^ and could 
be consistent with this inequality being satisfied as an 
equality- It also in very good agreement with the value of 
v — 0.69 obtained by Herbut^"*! in an expansion in powers 
away from the lower critical dimension di = 1- 

It has been suggestecl2Si^2i2S that the inequality ly > 
2/d can be violated and that in fact ly can be less than 
2/d- At the center of this debate is the way the average 
over the disorder is performed. If the disorder is cho- 
sen from a random distribution without any constraints, 
called the canonical ensemble of disorder, one might ask if 
equivalent results are obtained if constraints are imposed 
on the random distribution, the so-called micro-canonical 
ensemble of disorder. For instance, for the model consid- 
ered here one could constrain the random chemical poten- 
tial to have exactly the same average for each generated 
sample. The proof^^ of the "quantum" Harris criterion 
relies on the physically more relevant canonical ensemble 
of disorder being used. Subsequent work^Mi42iM44 have 
shown that even though the two ensembles in principle 
yield the same results in the thermodynamic limit, their 
finite-size properties can in some cases be different. All 
our results have been obtained using the canonical en- 
semble of disorder, in that at each site a potential was 
drawn from the uniform distribution [0, 1]. Hence, for a 
given sample the average chemical potential is not exactly 
1/2. In light of the above discussion we have therefore 
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performed additional simulations directly at Kc but this 
time imposing the constraint that the average chemical 
potential must be exactly 1/2 for every disorder realiza- 
tion. Our results for this micro-canonical ensemble of 
disorder are also shown in Fig|7|and are indistinguishable 
from our results obtained with the canonical ensemble. 
One should note that the uniform disorder distribution 
we employ is likely less sensitive to the difference between 
the canonical and micro-canonical ensemble of disorder 
than a bimodal distribution. Hence we conclude that our 
procedure for calculating v is valid. 

We finally discuss our results for the correlation func- 
tions which are shown in Fig. |S1 for a range of different 
system sizes calculated at K^- To determine the power 
law decay of the correlation functions we used lattices of 
size X aL^ with a = 1/8 and z = 3. All the temporal 
correlation functions can be fitted with the same expo- 
nent Hr = 1.06(7). We note that this value for yr would 
appear to satisfy the equality j/r < 1 as an equality. Fol- 
lowing section llil Al we assume that yr = {d + z — 2 + r])/z 
and using the above determined value for z = 3 we then 
infer 

ry^-1. (22) 

This would then satisfy Eq. (|15|l as an equality. Our 
results for the temporal correlation functions have been 
determined out to large lattice sizes Lj- = 512 and appear 
quite stable. We are therefore relatively confident that 
these results are trustworthy. 

Our results for the spatial correlation functions 
C(r,r — 0) appear less clear. The value we obtain 
for yr — 4.4(4) is clearly not consistent with our other 
estimates for z, since the ratio yr/yr = z then yields 
z = 4.15. If one assumes that this is the correct value 
for z it again follows that 77 ~ —1. We strongly suspect 
that this value of yr is due to the relatively small spatial 
extent {L < 16) of the lattices used. For such small lat- 
tice sizes the correlation functions have likely not reached 
their asymptotic behavior. We have varied the strength 
of the disorder, and found the same critical behaviour 
at A = 0.6 and A — 1.0. One might attempt to reach 



larger lattice sizes for the spatial correlation functions by 
decreasing the aspect ratio dramatically. We have at- 
tempted this by performing calculations for a lattice of 
size 20^ X 40 at Kc (also shown in Fig. correspond- 
ing to an aspect ratio of a = 1/200. In this case the 
spatial correlation functions show clear evidence for an 
exponential decay probably due to a dimensional cross- 
over induced by the extremely small aspect ratio. In fact, 
it would seem to be more reasonable to increase a to a 
more optimal aspect ratio determined by the requirement 
that {£,/LY ~ ^.r/Lr. Due to the large value of z, this 
has proven impossible for the present study. 

V. CONCLUSION 

In the present paper we have shown strong numerical 
evidence that the dynamical critical exponent z at the 
BG-SF critical point is equal to the spatial dimension 
for the three-dimensional site-disordered Bose-Hubbard 
model. The relation z = d therefore continues to hold 
in d = 3. Results for the scaling of the stiffness versus 
K , as well as at = 0.190(1) versus L^, are consistent 
with z = d = 3. The compressibility is almost constant 
and independent of system size through the critical point 
consistent with this value for z. In addition we have ob- 
tained values for the critical exponents v — 0.70(12) and 
?7 ~ — 1, with the cautionary note that a reliable deter- 
mination of "q is impeded by the very small spatial extent 
of the lattices available. In light of the recent work by 
Bernardet et alm^, it would have been very interesting 
to analyze each disorder realization seperately, in order 
to determine a characteristic Ki for each sample. The 
scaling analysis should then be redone using Ki. Unfor- 
tunately, we were not able to perform such an analysis 
for this work. 
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